---
title: "Transformer une sortie `R` (moche) en table de résultats (belle) avec `modelsummary()`"
date: 02-03-2026
categories: ["R", "Modèle linéaire", "Fonction"]
tags: ["Technique", "R"]
description: "Dans un récent post, je proposais des fonctions basiques pour transformer une sortie `lm()` en table de résultats. D'autres ont déjà fait le travail avant moi, et bien mieux. C'est donc une bonne raison pour présenter le package `modelsummary`."
image: img_display.png
image-alt: "Table de régression formatée finale"
lang: fr
body-classes: indent-paragraphs
format:
html:
code-fold: true
code-summary: "Voir le code"
code-overflow: wrap
code-tools: true
code-block-bg: true
execute:
warning: false
error: false
---
<hr style="border:1px solid blue" />
# Point de départ
```{r, echo=FALSE}
#### Charger les packages ####
library(tidyverse);
library(modelsummary)
# Charger les données
df <- palmerpenguins::penguins %>%
mutate(
species = case_when(
species == "Adelie" ~ "Adélie",
species == "Chinstrap" ~ "Jugulaire",
species == "Gentoo" ~ "Papou"),
sex = case_when(
sex == "female" ~ "Femelle",
sex == "male" ~ "Mâle")) %>% # Traduire les données
rename(
masse = body_mass_g,
espece = species,
sexe = sex) # Renommer les colonnes
```
Dans [un récent billet](https://charlymarie.github.io/posts/03%20-%20lm_summary()/lm_summary().html){target='_blank'}, je montrais comment automatiser l'extraction de coefficients d'un modèle `lm()`. Je trouve que c'est toujours bien de travailler *à la main* : manipuler les données, comprendre comment elles sont construites, extraites, comment la table est organisée puis rendue, etc. Cela permet de comprendre ce qui se passe sous le capôt de la machine, et de gagner en autonomie si un package n'existe pas déjà pour faire ce que l'on souhaite faire.
Mais, si l'on est honnête, mes fonctions étaient fragiles et la façon de faire peu élégante.
Aussi, ces fonctions étaient bornées à `lm()`. Un `glm()`, `ordinal()` [@ordinal], ou autre `glm.nb()` [@MASS] auraient renvoyé une erreur, et il aurait fallu construire d'autres fonctions pour en extraire les résultats. Il se trouve que d'autres ont déjà fait ce travail de construction, notamment pour le package `modelsummary` [@modelsummary] !
Je repars du même jeu de données `penguins` [@palmerpenguins] et du même modèle linéaire :
$$
\text{masse}_i
= \alpha
+ \beta_1\,X_{1i}
+ \beta_2\,X_{2i}
+ \beta_3\,X_{3i}
+ \varepsilon_i
$$ {#eq-un}
où
$$
X_{1}=\mathbb{1}(\text{espèce}_i=\text{Jugulaire});\quad
X_{2}=\mathbb{1}(\text{espèce}_i=\text{Papou});\quad
X_{3}=\mathbb{1}(\text{sexe}_i=\text{Mâle}).
$$
A nouveau, la sortie en base `R` n'est pas très élégante.
```{r, output=axis, echo=FALSE}
#| label: tbl-regression-moche
#| tbl-cap: "Régression de la masse des manchots en grammes, sur l'espèce et le sexe, selon l'@eq-un. Sortie base `R`."
lm_1 <- lm(masse ~ espece + sexe, df)
summary(lm_1)
```
# La magie de `modelsummary::modelsummary()`
`modelsummary()` va être très pratique pour un très grand nombre de raisons :
- Supporte énormément de modèles, inférentiels ou bayésiens, des plus obscurs aux plus importants [comme `fixest`, @berge2026].
- Directement combinable avec des erreurs personnalisées, du type robustes, clusterisées, etc, via `sandwich` [@sandwich].
- Largement personnalisable, comme on va le voir.
- Les tables sont élégantes.
Dans sa version la plus basique, la fonction renvoie une table de régression comme on en voit tant.
```{r}
#| label: tbl-regression-modelsummary1
#| tbl-cap: "Régression de la masse des manchots en grammes, sur l'espèce et le sexe, selon l'@eq-un. Sortie basique de `modelsummary()`."
modelsummary(lm_1)
```
En soit, c'est tout ce dont on a besoin : à partir de l'estimation et de son erreur, on peut approximer *t*, *p*, et l'intervalle de confiance à 95 %. Mais, l'on pourrait tout de même souhaiter afficher les véritables statistiques plutôt que de devoir les approximer.
J'ajoute donc les statistiques suivantes, alignées dans cet ordre dans la @tbl-regression-modelsummary2 :
- Estimation.
- Intervalles de confiance à 95 %.
- Erreur standard.
- Statistique *t*.
- Statistique *p*.
```{r}
#| label: tbl-regression-modelsummary2
#| tbl-cap: "Régression de la masse des manchots en grammes, sur l'espèce et le sexe, selon l'@eq-un. Sortie `modelsummary()` avec statistiques."
modelsummary(lm_1,
estimate = c("Estimation <br> [IC 95 %]" = "{estimate}<br>[{conf.low}; {conf.high}]"), # Ajouter l'intervalle de confiance à l'estimate
statistic = c(
"Erreur standard" = "std.error",
"Statistique t" = "statistic",
"p-valeur" = "p.value")) # Ajouter la statistique t et la p-valeur
```
Ici, ça se complique : les statistiques s'empilent, la table perd en lisibilité. On va devoir la faire pivoter, pour représenter les mêmes statistiques en colonne plutôt qu'en ligne.
```{r}
#| label: tbl-regression-modelsummary3
#| tbl-cap: "Régression de la masse des manchots en grammes, sur l'espèce et le sexe, selon l'@eq-un. Sortie `modelsummary()` avec statistiques pivotées."
modelsummary(lm_1,
estimate = c("Estimation <br> [IC 95 %]" = "{estimate}<br>[{conf.low}; {conf.high}]"),
statistic = c(
"Erreur standard" = "std.error",
"Statistique t" = "statistic",
"p-valeur" = "p.value"),
shape = term ~ statistic) # Pivoter la table
```
Cette fois, la table commence à être exploitable. Je finalise l'aspect cosmétique, en :
1. Renommant les coefficients.
2. Retirant l'ensemble des indices d'ajustement, en renvoyant uniquement le nombre d'observations et le *R²* ajusté en une note de bas de tableau.
3. Retirant le numéro de modèle (1).
4. Mettant les titres de colonnes en *gras* et les statistiques en *italique*.
5. Formattant les statistiques à deux décimales, et la *p*-valeur selon les normes APA.
6. Ajoutant le nombre de degrés de liberté à la colonne *t*
```{r}
#| label: tbl-regression-modelsummary4
#| tbl-cap: "Régression de la masse des manchots en grammes, sur l'espèce et le sexe, selon l'@eq-un. Sortie `modelsummary()` avec statistiques cosmétiques."
g <- broom::glance(lm_1) # Stocker les indices d'ajustement
ddl <- unique(get_estimates(lm_1)$df.error) # Stocker les degrés de liberté
modelsummary(
shape = term ~ statistic,
estimate = c(
"<strong>Estimation<br>[IC 95 %]</strong>" = "{estimate}<br>[{conf.low}; {conf.high}]"), # <strong> pour mettre en gras et <em> pour mettre en italique
statistic = setNames(
c("std.error", "statistic", "p.value"),
c(
"<strong>Erreur standard</strong>",
paste0("<strong>Statistique <em>t</em> (", ddl, ")</strong>"), # On vient ici ajouter les degrés de liberté à la colonne
"<strong><em>p</em>-valeur</strong>")),
models = list(" " = lm_1), # Supprimer le (1)
gof_omit = ".*", # Supprimer les indices d'ajustement
notes = (
paste0(
"<em>Note</em> : <em>n</em> = ",
scales::number(g$nobs),
" observations. <em>R²</em> ajusté = ",
sprintf("%.2f", g$adj.r.squared),
".")), # Ajouter les indices d'ajustement en note de bas de tableau
coef_rename = c(
"especeJugulaire" = "Espèce : Jugulaire",
"especePapou" = "Espèce : Papou",
"sexeMâle" = "Sexe : Mâle"), # Renommer les coefficients
fmt = fmt_statistic(estimate = 2, std.error = 2, statistic = 2, conf.low = 2, conf.high = 2, p.value = scales::label_pvalue())) # Formatter les décimales
```
# Construire une fonction personnalisée et simple, regroupant tous les besoins classiques en psychologie : `fr_modelsummary_lm()`
Il ne manque finalement qu'une chose à cette table : d'inclure un coefficient standardisé $\beta$. Cela demande un peu plus de travail car, de base, `modelsummary()` ne renvoie pas de coefficient standardisé. On va donc devoir (à nouveau) créer une table personnalisée avant de pouvoir l'envoyer à la fonction.
Plutôt que de devoir à chaque fois ré-écrire ces scripts, et potentiellement devoir les mettre à jour en cas de changement, je construits une fonction pour (1) extraire une table personnalisée et (2) formater cette table personnalisée. Je m'appuie encore une fois sur le package `parameters` [@parameters] pour simplifier l'importation d'une matrice de variance-covariance facilement personnalisable^[Pour les personnes un peu plus geek, il serait aussi possible de [définir une fonction personnalisée et directement utilisable en interne de `modelsummary()`](https://modelsummary.com/vignettes/modelsummary_extension.html#modifying-information-tidy_custom-and-glance_custom){target='_blank'}, de type `tidy_custom`.].
Comparées au post précédent, ces fonctions sont directement intégrées l'une à l'autre et permettent d'appeler une seule fonction personnalisée `fr_modelsummary_lm()` (ou son équivalent anglophone `en_modelsummary_lm()`, que je vous mets plus bas). Elles acceptent également un argument `vcov`. Finalement, elles seraient facilement personnalisables et adaptables à d'autres modèles avec quelques ajustements.
```{r}
# Fonction fr
fr_modelsummary_lm <- function(
model,
vcov = NULL,
coef_rename = NULL,
gof_omit = ".*",
...) {
# Vérifier que le modèle est bien de type lm()
if (!inherits(model, "lm")) {
stop("Le modèle fournit doit être un objet de classe lm. Ce n'est pas le cas ici.",
call. = FALSE)
}
# Première fonction : Extraire une table
custom_modelsummary_lm <- function(model, standardize = "refit", vcov = NULL, ...) {
est <- modelsummary::get_estimates(model, vcov = vcov) # Extraire les coefficients pour la table de base
est_std <- parameters::model_parameters(model, standardize = standardize, vcov = vcov) %>% # Extraire un coefficient standardisé
dplyr::select(Parameter, Coefficient, CI_low, CI_high) %>% # Sélectionner les colonnes standardisées
dplyr::rename(Std_Coefficient = Coefficient) # Renommer la colonne qui nous intéresse
est$estimate_std <- paste0(
sprintf("%.2f", est_std$Std_Coefficient), "<br>",
"[", sprintf("%.2f", est_std$CI_low), "; ", sprintf("%.2f", est_std$CI_high), "]") # Formatter les statistiques à deux décimales
out <- list(
tidy = est,
glance = broom::glance(model)) # Stocker la table et les indices de fit
class(out) <- "modelsummary_list" # Définir la classe de l'objet de telle façon à ce que modelsummary() l'accepte
out
}
# Seconde fonction : Formatter la table
est <- custom_modelsummary_lm(model, vcov = vcov) # Appliquer la première fonction au modèle
g <- broom::glance(model) # Stocker les indices d'ajustement
ddl <- unique(modelsummary::get_estimates(model, vcov = vcov)$df.error) # Stocker les degrés de liberté
vcov_note <- if (!is.null(vcov)) paste0(" Erreurs standards : <em>", vcov, "</em>.") else ""
notes_txt <- paste0(
"<em>Note</em> : <em>n</em> = ",
scales::number(g$nobs),
". <em>R²</em> ajusté = ",
sprintf("%.2f", g$adj.r.squared),
".",
vcov_note)
modelsummary(
coef_rename = coef_rename,
gof_omit = gof_omit,
shape = term ~ statistic,
estimate = c(
"<strong>Estimation<br>[IC 95 %]</strong>" = "{estimate}<br>[{conf.low}; {conf.high}]"),
statistic = setNames(
c("std.error", "statistic", "p.value", "estimate_std"), # Ajouter la colonne standardisée
c(
"<strong>Erreur standard</strong>",
paste0("<strong>Statistique <em>t</em> (", ddl, ")</strong>"),
"<strong><em>p</em>-valeur</strong>",
"<strong><em>β</em><br>[IC 95 %]</strong>")), # Renomme
models = list(" " = est),
notes = notes_txt,
fmt = fmt_statistic(estimate = 2, std.error = 2, statistic = 2, conf.low = 2, conf.high = 2, p.value = scales::label_pvalue()),
...)
}
```
Ces fonctions personnalisables s'appuient ultimement sur `modelsummary()` et tout argument supplémentaire lui sera donc passé. Un exemple en français via `fr_modelsummary_lm()`, en intégrant une matrice de variance-covariance robuste à l'hétéroscédasticité.
```{r}
#| label: tbl-regression-modelsummary5
#| tbl-cap: "Régression de la masse des manchots en grammes, sur l'espèce et le sexe, selon l'@eq-un. Sortie `fr_modelsummary_lm()`."
fr_modelsummary_lm(lm_1,
coef_rename = c(
"especeJugulaire" = "Espèce : Jugulaire",
"especePapou" = "Espèce : Papou",
"sexeMâle" = "Sexe : Mâle"),
vcov = "HC3")
```
# La fonction `en_modelsummary_lm()`
Pour celles et ceux qui pourraient en avoir besoin en anglais, la fonction traduite et le même exemple en appelant cette fois la fonction `en_modelsummary_lm()` :
```{r}
#| label: tbl-regression-modelsummary6
#| tbl-cap: "Régression de la masse des manchots en grammes, sur l'espèce et le sexe, selon l'@eq-un. Sortie `en_modelsummary_lm()`."
# English function
en_modelsummary_lm <- function(
model,
vcov = NULL,
coef_rename = NULL,
gof_omit = ".*",
...) {
# Check that the model is of type lm()
if (!inherits(model, "lm")) {
stop("The model must be an object of class lm. This is not the case here.",
call. = FALSE)
}
# First function: Extract a table
custom_modelsummary_lm <- function(model, standardize = "refit", vcov = NULL, ...) {
est <- modelsummary::get_estimates(model, vcov = vcov) # Extract coefficients for the base table
est_std <- parameters::model_parameters(model, standardize = standardize, vcov = vcov) %>% # Extract standardised coefficients
dplyr::select(Parameter, Coefficient, CI_low, CI_high) %>% # Select standardised columns
dplyr::rename(Std_Coefficient = Coefficient) # Rename the column we are interested in
est$estimate_std <- paste0(
sprintf("%.2f", est_std$Std_Coefficient), "<br>",
"[", sprintf("%.2f", est_std$CI_low), "; ", sprintf("%.2f", est_std$CI_high), "]") # Format statistics to two decimal places
out <- list(
tidy = est,
glance = broom::glance(model)) # Storing the table and fit indices
class(out) <- "modelsummary_list" # Define the class of the object so that modelsummary() accepts it
out
}
# Second function: Format the table
est <- custom_modelsummary_lm(model, vcov = vcov) # Apply the first function to the model
g <- broom::glance(model) # Stock fit indices
ddl <- unique(modelsummary::get_estimates(model, vcov = vcov)$df.error) # Stock degrees of freedom
vcov_note <- if (!is.null(vcov)) paste0(" Standard errors : <em>", vcov, "</em>.") else ""
notes_txt <- paste0(
"<em>Note</em> : <em>n</em> = ",
scales::number(g$nobs),
". Adjusted <em>R²</em> = ",
sprintf("%.2f", g$adj.r.squared),
".",
vcov_note)
modelsummary(
coef_rename = coef_rename,
gof_omit = gof_omit,
shape = term ~ statistic,
estimate = c(
"<strong>Estimate<br>[95% CI]</strong>" = "{estimate}<br>[{conf.low}; {conf.high}]"),
statistic = setNames(
c("std.error", "statistic", "p.value", "estimate_std"), # Add standardised column
c(
"<strong>Standard error</strong>",
paste0("<strong><em>t</em> statistic (", ddl, ")</strong>"),
"<strong><em>p</em>-value</strong>",
"<strong><em>β</em><br>[95% CI]</strong>")), # Rename column
models = list(" " = est),
notes = notes_txt,
fmt = fmt_statistic(estimate = 2, std.error = 2, statistic = 2, conf.low = 2, conf.high = 2, p.value = scales::label_pvalue()),
...)
}
en_modelsummary_lm(lm_1,
coef_rename = c(
"especeJugulaire" = "Species : Chinstrap",
"especePapou" = "Species : Gentoo",
"sexeMâle" = "Sex : Male"),
vcov = "HC3")
```
# Références